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Abstract 

Particle-number projection within the Lipkin-Nogami (LN) method is applied to the self- 
consistent quasiparticle random-phase approximation (SCQRPA), which is tested in an exactly 
solvable multi- level pairing model. The SCQRPA equations are numerically solved to find the 
energies of the ground and excited states at various numbers of doubly degenerate equidistant 
levels. The use of the LN method allows one to avoid the collapse of the BCS (QRPA) to obtain 
the energies of the ground and excited states as smooth functions of the interaction parameter G. 
The comparison between results given by different approximations such as the SCRPA, QRPA, 
LNQRPA, SCQRPA and LNSCQRPA is carried out. While the use of the LN method significantly 
improves the agreement with the exact results in the intermediate coupling region, we found that 
in the strong coupling region the SCQRPA results are closest to the exact ones. 
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I. INTRODUCTION 



The random-phase approximation (RPA), which includes correlations in the ground state, 
provides a simple theory of excited states of the nucleus. However, the RPA breaks down at 
a certain value G cr of interaction parameter G, where it yields imaginary eigenvalues. The 
reason is that the RPA equations, linear with respect to the X and Y amplitudes of the 
RPA excitation operator, are derived based on the quasi-boson approximation (QBA). The 
latter neglects the Pauli principle between fermion pairs and its validity is getting poor with 
increasing the interaction parameter G. The collapse of the RPA at the critical value G CT of 
G invalidates the use of the QBA. The RPA therefore needs to be extended to correct this 
deficiency, at least for finite systems such as nuclei. 

One of methods to restore the Pauli principle is to renormalize the conventional RPA 
to include the non-zero values of the commutator between the fermion-pair operators in 
the correlated ground state. These so-called ground-state correlations beyond RPA are ne- 



glected within the QBA. The interaction in this way is renormalized and the col 



apse of RPA 
f 

3[. However, 



is avoided. The resulting theory is called the renormalized RPA (RRPA) 
the test of the RRPA carried out within several exactly solvable models showed that the 
RRPA results are still fax from the exact solutions 0, 0, B|. 

Recently, a significant development in improving the RPA has been carried out within 
the self-consistent RPA (SCRPA) Based on the same concept of renormalizing 

the particle-particle (pp) RPA, the SCRPA made a step forward by including the screening 
factors, which are the expectation values of the products of two pairing operators in the 
correlated ground state. The SCRPA has been applied to the exactly solvable multi-level 
pairing model, where the energies of the ground state and first excited state in the sys- 
tem with N + 2 particles relative to the energy of the ground-state level in the iV-particle 
system are calculated and compared with the exact results. It has been found that the 
agreement with the exact solutions is good only in the weak coupling region, where the 
pairing-interaction parameter G is smaller than the critical values G cr . In the strong cou- 
pling region {G » G cr ), the agreement between the SCRPA and exact results becomes 
poor 4, [5( • In this region a quasiparticle representation should be used in place of the pp 
one, as has been pointed out in Ref. [z|. As a matter of fact, an extended version of the 
SCRPA in the superfluid region has been proposed and is called the self-consistent quasipar- 
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ticle RPA (SCQRPA), which was applied for the first time to the seniority model in Ref. 
and a two-level pairing model in Ref. [sj]. However, the SCQRPA also collapses at G — G CI . 
It is therefore highly desirable to develop a SCQRPA that works at all values of G and also 
in more realistic cases, e.g. multi-level models. The aim of the present work is to construct 
such an approach. Obviously, the collapse of the SCQRPA at G = G CT , which is the same 
as that of the non-trivial solution for the pairing gap within the Bardeen-Cooper-Schrieffer 
theory (BCS), can be removed by performing the particle-number projection (PNP). The 



Lipkin-Nogami method [10|, which is an approximated PNP before variation, will be 
used in such extension of the SCQRPA in the present paper because of its simplicity. This 
approach shall be applied to a multi-level pairing model, the so-called Richardson model 
12j, which is an exactly solvable model extensively employed in literature to test approxi- 
mations of many-body problems. 

The paper is organized as follows. Section [III presents a brief outline of the SCQRPA 
theory that includes the PNP within the LN method. The results of numerical calculations 
are analyzed and discussed in sec. IHII Conclusions are drawn in the last section. 

II. FORMALISM 

A. Model Hamiltonian 



The Richardson model (also called the multi-leve 



)airing 



model, picket-fence model or 



121 ]. It consists of Q doubly-fold 



ladder model, etc) was described in detail in Refs. {4, 0, [f], 
equidistant levels interacting via a pairing force with a constant parameter G. The model 
Hamiltonian is given as 

n n 

H = $> ~ W - G E p J p r > C 1 ) 

j=i i,i'=i 
where €j are the single-particle energies on the j-shells. The particle-number operator Nj 
and pairing operators Pj , Pj on the j-th orbital (with unit shell degeneracy j + 1/2 = 1) 
are defined as 

Nj = a'jCLj + a[_ j a_ j , (2) 

v.- Pi = (p}) ] - (3) 
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These operators fulfill the following exact commutation relations 

[Pi , P}] = 8 jr (l-N j ) , (4) 
[Nj , Pj,} = 28 jr Pj l , [Nj , P r ] = -28 jf P r . (5) 

By using the Bogoliubov transformation from particle operators a, and aj to quasiparticle 
ones a] and aj 

aj = UjCtj + VjCt-j , a^j = Uja-j — VjOtj , (6) 
he pairing Hamiltonian in Eq. (CQ) is transformed into the quasiparticle Hamiltonian as 



13, 



i i ii' ii' 

+ ^ (^UJ, + A f Aj) + ^ ^/A' ; AV , (7) 

ii' ii' 

where A/} is the quasiparticle-number operator, while Aj and Aj are the creation and de- 
struction operators of a pair of time-conjugated quasiparticles: 

Afj = OjOCj + aljOc-j , (8) 

The commutation relations between operators A/}, and Aj are similar to those for particle 
operators in Eqs. (jl]) and ([5]), namely 

[A- , = foil - Ms) > ( 10 ) 
[Aj , 4] = , [Nj , A jf ] = • (11) 

The coefficients a, bj, Cj, djj>, gj(f), hjj>, qjj> in Eq. (j7|) are given in terms of the coefficients 
Uj, Vj of the Bogoliubov transformation, and the single particle energies €j as (see, e.g. Ref. 
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3 \ 3 J 3 



(12) 



h 3 = ( £ J - A )( M i - V ]) + ^GUjVj U 3' V 3' + Gv j , ( 13 ) 

r 

= 2(ej - X)UjVj - G{u) - u?) ^ W , - 2GVJ , (14) 

i' 
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djji = —G(u?u?, + v?Vj,) = dfj , (15) 
9j(f) = GujVj(Uj, - v 2 jt ) , (16) 

K 1 = jitfVj* + v2 j U )') = k fj > ( 17 ) 

Qjj' = -GujVjUfVj, = q fj . (18) 

The single-particle energies are given as ej = je, where j = 1, . . . , Q, and e is the level 
distance chosen to be equal to 1 MeV in the present work. The chemical potential A and 
the coefficients Uj and Vj are determined by solving the gap equations discussed in the next 
section. 



B. Gap and number equations 

1. Renormalized BCS 

It is well known that the Pauli principle between the quasiparticle-pair operators Aj and 
A]j, is neglected within the conventional BCS, which assumes that (BCS|A/}|BCS) = within 
the BCS ground state |BCS). A simple way to restore the Pauli principle is to introduce a 
new ground state |0) in which the correlations among quasiparticles lead to nonzero values 
of the quasiparticle occupation numbers so that the contribution of the A/j-term at the right- 
hand side (rhs) of Eq. fflOl) is preserved. By doing so, the BCS equations are renormalized 
and the resulting theory is called the renormalized BCS (RBCS) [15|. Within the RBCS the 
commutator between the quasiparticle-pair operators are defined as 

(0|[^,4P>=MZ>i>, (19) 

with 

V s = l-N Si (V,) = l-2n 3 , (20) 
where rij is the quasiparticle number in the correlated ground state |0) 

^. = 1(01^10)^0. (21) 

Taking into account Eq. (1191) and performing a constrained variational calculation to min- 
imize the Hamiltonian H = H' — \N, where iV = . Nj is the particle-number operator, 



5 



the RBCS equations for the pairing gap A and particle number N have been derived as [15J 

A = Gj>' N = 2^ Pj , (22) 



where 



1 



Tj = UjVjpj) , Pi = vj(Vj) + -(1 - (T>j)) , (23) 



^' = \/(ej-Gv]-X) 2 + A 2 . (25) 

The renormalization factors (T>j), called the ground-state correlation factors, are obtained by 
solving the SCQRPA equations discussed later in this paper (See Sec. Ill D 2\i . The internal 
energy of the system within the RBCS ground state (the RBCS ground-state energy) is 
given as 

E RBCS = 2 £ (c . _ x)p . -^-GJ2pI (26) 

3 3 

By setting (T>j) = 1, the RBCS equations go back to the well-known BCS ones. 



2. BCS with SCQRPA correlations 

In the minimization procedure, which leads to the equation (See, e.g. Ref. [if] ]) 

(0|[F,^]|0) = 0, (27) 

the RBCS ignores the expectation values (A],Aj) = (0\A^Aj\0) and (AfAj) = (O^^lO) 
of the products of pair operators in the correlated quasiparticle ground state |0). By retaining 
these screening factors in calculating the left-hand side (lhs) of Eq. (127]) . we derive from Eq. 
(127]) an equation for the level-dependent pairing gap in the form 

A? = C- J \.^ } , (28) 

with the single-particle energies ej in the expressions for Uj and Vj in Eq. (124]) being renor- 



malized to e'j as 



; ~ ^ -r E(4 - 4) (<44> + (A'A')) • (29) 



G 

e,- = e,- + 
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We call Eq. (I28I) the BCS gap equation with SCQRPA correlations, and use the abbreviation 
BCS1 to denote this approach, having in mind that it includes the screening factors (A-A|,) 
and {AjAj') in the renormalized single-particle energies given by Eq. (1291) . These screening 
factors are found by solving Eqs. f[2"5j) and f[2"9"j) selfconsistently with the SCQRPA ones to 
be discussed later in Sec. IHD[ where the explicit expressions of the screening factors are 
given in terms of the SCQRPA forward- and backward-going [X and y) amplitudes. The 
limit case of Eqs. fl28|) and fl29|) for a degenerate two-level model is studied in Ref . js] . 

The rhs of Eq. (|28|) contains the expectation values (VjVji), whose exact treatment is 
not possible as it involves an infinite series in terms of the products of A^AjA^Af {9], 
or an infinite boson expansion [17], which again needs to be truncated at a certain order. 



In Ref. (9] this series is truncated at the first order, while the consideration in Ref. 17] 
is limited up to the four-boson terms. Such expansion is based on the method of treating 
the single-particle (quasiparticle) density used by Rowe in Ref. Q] or a mapping employed 
in Ref. [3j. In the numerical calculations within the present paper we treat these terms 
approximately as follows. By noticing that the expectation values (T>jT>j/) are present in in 
the ratios (VjVj/) / (Vj) or, more general, (VjVj/) / a/ (Vj) (Dp), and that 

PPy) = (Vj) (V f ) + 5Af jf , with 8Njj> = <JW r ) - (Afj) (Nj>) , (30) 

we rewrite these ratios as 



(VMVj,) + , " . (31) 

v #W 

The numerator 8J\fjj> of the last term at the rhs of Eq. (!3T|) can be estimated by using the 
mean-field contraction as 

5N jr - 2^(1 - //,) = StfiSNj) 2 , (32) 

where (SAfj) 2 = (A/J) — (ftfj) 2 = 2n 3 -(l — rij) is the quasiparticle-number fluctuation on the 
j-th orbital. This quantity is much smaller than 1, while the denominator y {Dj)(T>f), 
which is also the first term at the rhs of Eq. (131j) . is comparable with 1 as the ground-state 
correlations factors (Vj) are not much smaller than 1. Therefore the last term at the rhs of 
Eq. (j3Tj) can be safely neglected so that 

(vp r ) ^ 



(V 3 )(V f ) . (33) 



Consequently, the ratio (T>jT>j>) / (T>j) in the sum over j' at the rhs of Eq. (|28|) can be simply 
approximated with (T>ji) 1 . In this case Eq. fl28|) takes the same level-independent form as 
that of Eq. (1221) for the RBCS gap except that the single-particle energies in uy and vj> are 
now given by Eq. (12~§1) . In the rest of the paper, such level- independent approximation for 
the pairing gap is assumed, whose numerical accuracy is checked in the Appendix [A3 



C. Lipkin-Nogami method with SCQRPA correlations 

The main drawback of the BCS is that its wave function is not an eigenstate of the 
particle- number operator N. The BCS, therefore, suffers from an inaccuracy caused by 
the particle-number fluctuations. The collapse of the BCS at a critical value G cr of the 



pairing parameter G, below which it has only a trivial solution wit 



i zero pairing gap, is 



intimately related to the particle-number fluctuations within BCS This defect is cured 
by projecting out the component of the wave-function that corresponds to the right number 
of particles. The Lipkin-Nogami (LN) method is an approximated PNP, which has been 



shown to be simple and yet efficient in many realistic calculations (See Ref. [18j for a 
recent detailed clarification of the use of the LN method). This method, discussed in detail 

nn 

Refs. [10|, HJJ, is a PNP before variation based on the BCS wave function, therefore the 
Pauli principle between the quasiparticle-pair operators Eq. (imp is still neglected within the 
original version of this method. In the present work, to restore the Pauli principle we propose 
a renormalization of the LN method, which we refer to as the renormalized LN (RLN) 
method or LN method with SQRPA correlations (LN1) when they are based on the RBCS or 
BCS1, respectively. Similar to the BCS1 (RBCS), the LN1 (RLN) includes the quasiparticle 
correlations in the correlated ground state |0), and the LN1 (RLN) equations are obtained 
by carrying out the variational calculation to minimize Hamiltonian H = H' — XN — X2N 2 . 
The LN1 equations obtained in this way have the form 

£ = G2>, N = 2J2pj, (34) 



1 In Refs. 0,0] the factorization (NjNf) ~ {Nj){Nj>) (j — p, h) was straightforwardly used to close the 
SCRPA equations because (Nh)(Nh'), whose value in the Hartree-Fock (HF) limit is 4, is much larger 
than the particle-number fluctuation (SN^) 2 = 2/^(1 — fh)- This is no longer the case for quasiparticle 
numbers, where (SNj) 2 are of the same order with (Afj) 2 ■ 
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€j = e 3 + (4A 2 - G)v) , A = Ai + 2A 2 (A + 1) 



(35) 



where 



rj = ujVjiVi) , pj = v*(Vj) + -(1 - (Vj)) 



< ' 1 V 



(6, - A)' + 



(36) 
(37) 



The coefficient A 2 has the following form 19] 



G E,(i - E," fro - E,(i - p,) 2 p1 



E,Pi(i-Pi) -E 7 (i-p,) 2 ^ 



(38) 



which becomes the expression given in the original paper [ll[ of the LN method in the limit 
of (T>j) = 1 and e'j = €j. The internal energy obtained within the LN1 ground state (the 
LN1 ground-state energy) is given as 

2 E ^ - X ^ ~ ^ " G E P) - A 2 AiV2 > ( 39 ) 



LNl 



where the expression for the particle-number fluctuation AN 2 in terms of uj, Vj and n 



(1 — (Dj))/2 has been derived in Ref. 14j . The LNl equations becomes the RLN equations 
by replacing the renormalized single-particle energies e'- defined in Eq. (1291) with ej. The 
RLN equations return to the BCS ones in the limit case, when A 2 = and (T>j) = 1. 



D. SCQRPA equations 



1. QRPA 



The QRPA excited state \u) is constructed by acting the QRPA operator Q\ 



Qt = YJWA ~ Y i A i) > Q» = (Q 



vl 1 



(40) 



on the QRPA ground state |0) as 

W)=Ql\o), 

where |0) is defined as the vacuum for the operator (14*01) . i.e. 



(41) 



Q„|0) = . 



(42) 
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The QBA assumes the following relation 



(0|[^,4]|0> = % 



(43) 



Within the QBA the QRPA amplitude Xj and obey the well-known normalization 
(orthogonality) conditions 



(44) 



(45) 



to guarantee that the QRPA operators (1401) are bosons, i.e. 

(0| \Qu,Ql] \0) = 5 m > . 



By linearizing the equation of motion with respect to Hamiltonian (jxj) and operators ( 1401 ). 
the set of linear QRPA equations is derived and presented in the matrix form as follow 




-Yf 



(46) 



where the QRPA submatrices are given as 



2(b j + 2q jf )6 jf + d 



33 i 



B 33' ~ 2 ( 1 ~ $33') h 33' ■ 



(47) 
(48) 



and the eigenvalues u v = S v — £q are the energies E v of the excited states relative to that of 
the ground-state level, £ . The QRPA ground-state energy is given as the sum of the BCS 



ground-state energy E^ s and the QRPA correlation energy as follows 



2Q| 



£QRPA = E _ 



BCS 



(49) 



2. Renormalized QRPA 



To restore the Pauli principle, the QRPA is renormalized based on Eq. (I19p instead of 
the QBA (1431) . The RQRPA operators are introduced as [20J] 



(50) 
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which are bosons within the quasiparticle correlated ground state |0), i.e. 



Qu, Ql 



|0) = <W , 



(51) 



if the X? and y*f amplitudes satisfy the same orthogonality conditions ffl4|) . namely 

£ (*T< " = <W ■ (52) 

The RQRPA submatrices are given as 



Ajjf — 2 



6 . + 2? .. / + 2 ^ ? ..,,^i_i^dj 



(VP?) 



(53) 
(54) 



The ground-state correlation factor (T>j) has been derived as a function of the backward- 
going amplitudes y^ (see e.g. Refs. |3|, |20|]) as 

1 



1 + £„ (y/) 



2 ' 



(55) 



whose values are found by consistently solving Eq. fl55l) with the RQRPA equations under 
the orthogonality condition (jUJ) for X^ and y 1 - amplitudes. In the limit of (Vj) = 1, one 
recovers from Eqs. ([531), ($M§ the QRPA matrices (03) and ( HHl) . 



5. SCQRPA and Lipkin-Nogami SCQRPA 



The only difference between the SCQRPA and the RQRPA is that, similarly to the 
SCRPA 0, fl, 0, the SCQRPA includes the screening factors, which are the expectation 
values of the pair operators (Aj,Aj) and (Aj'Aj) over the correlated quasiparticle ground 
state |0). The SCQRPA operators are defined in the same way as that for the RQRPA ones 
so is the correlated ground state. Therefore we use for it the same notation |0) having in 
mind the above-mentioned difference due to screening factors. 

The SCQRPA submatrices are obtained in the following form 

1 



.4 



j.i 



<Pj) 



£ (l >r < A" A-> - 2 £ h »" (AfAj) 



33' 



(56) 
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+d n' r^rri^ + Hi,' 



hjf + -^-y ( E djfiAfAj) + 2 E h 33"(4" A ^ 



3" 3 



+ 2hjj , jPPf) }4*r) , (57 ) 

where the screening factors (At Ay) and (AjAf) are given in terms of the amplitudes Xj 
and y V A as 



<^U,v) = (0|^|0> = ^(P^P , (58) 



(A- A-') = (o|AA"|o) = y/CDMVf) Y, x Pr ■ ( 59 ) 

V 

The rhs of Eqs. (I58|) and are obtained by using the inverted transformation of Eq. (1501) . 
namely 



4 = v^EC^+^O ' ( 6 °) 

and Eq. flBll . 

For the internal (ground-state) energy, the relation ( )49|) no longer holds due to the pres- 
ence of the ground-state correlation factors (T>j) in the SCQRPA equations. Therefore, the 
SCQRPA ground-state energy is calculated directly as the expectation value of the Hamil- 
tonian ([7j) in the correlated quasiparticle ground state, namely 

^scQRPA = m - 0) = a + J2 6,(1 - (V 3 )) + d jf (A]A f ) 

3 33' 

+ h 33> ((44') + (A'A)) + E - D i) - > ( 61 ) 



In the numerical calculations in the present paper the exact ratios (DjDji) / <J (T>j) (T>ji) 
in the RQRPA and SCQRPA submatrices fl53l) . (J5lj) . fl56|) . and fl57|) are calculated within 
the approximation ( 1331) . whose accuracy within the SCQRPA is numerically tested in the 
Appendix [Al 

Concerning the SCQRPA ground-state energy, by using Eq. fl30|) and relation fl32|) . the 
last term at the rhs of Eq. (16T|) can be approximated as 

E %•'((! - ^) (i - %)) ^ E^ ^ - to)) ( x - to)) + E^^) 2 

ii' ii' i 
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GA 



(1 _ (p.)) (i _ (p.,)) i 



(62) 



J J' 3 3 3 3 

The set of Eq. (1241) (for Uj and Uj) with the renormalized single-particle energies e'j ( 1291) 
replacing e,,-, Eq. (j4"6"l) with submatrices (156|) . (137)1 . and Eq. (1521) (for the amplitudes AJ, 

and energies u v ), together with Eq. (l55l) (for the ground-state correlation factors (T>j)) 
forms a set of coupled non-linear equations for Uj, Vj, X^, y?, u u , and (T>j). This set is 
solved by iteration in the present paper to ensure the self-consistency with the SCQRPA. 
Neglecting the screening factors ([SHD and (1511 the SCQRPA is reduced to the RQRPA, and 
the SCQRPA correlated ground state |0) becomes the RQRPA ground state. 

The Lipkin-Nogami SCQRPA (LNSCQRPA) equations have the same form as that of the 
SCQRPA ones given in Eqs. (156]) and (l5"Tj) . but the chemical potential and coefficients of 
the Bogoliubov transformation are determined by solving the LN1 gap equations (l34"j) . (|35|) 
instead of the BCS ones. 



III. ANALYSIS OF NUMERICAL CALCULATIONS 

We carried out the calculations of the ground-state energy, Eg_ s , and energies of excited 
states, w v = £ v — S o , in the quasiparticle representation using the BCS, QRPA, SCQRPA 
as well as their renormalized and PNP versions, namely the RBCS, BCS1, LN, RLN, LN1, 
LNQRPA, and LNSCQRPA, at several values of particle number N. The detailed discussion 
is given for the case with = 10. In the end of the discussion we report a comparison 
between results obtained for = 4, 6, 8, and 10 to see the systematic with increasing N. 



A. Pairing gap 

Shown in Fig. [T]are the pairing gaps obtained within the BCS, RBCS, BCS1, LN, RLN, 
and LN1 as functions of the pairing-interaction parameter G for = 10. Similarly to the 



two-level case 20J], the BCS has only a trivial solution Abcs = at G < G^ cs = 0.34 
MeV, while at G > G* cs the gap A BC s increases with G. Within the BCS1 (RBCS) the 
ground-state correlation factor (Dj) is always smaller than 1 (at G / 0). This shifts up 
the value of the critical point G CI to G^ r BCS ~ 0.38 MeV, and G BCS1 ~ 0.47 MeV so that 
G BCS < G^ r BCS < G B r csl . The PNP within the LN method completely smears out the BCS 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 

G (MeV) 

FIG. 1: (Color on line) Pairing gaps A as functions of G for iV=i0. The dotted, thin and thick 
dash-dotted denote the BCS, RBCS, and BCS1 results, respectively while the dashed, thin and 
thick dash-double-dotted lines represent the LN, RLN, and LN1 results, respectively. 

and BCS1 (RBCS) critical points to produce the pairing gap Aln as a smooth function of 
G, which increases with G starting from its zero value at G = 0. It is worth noticing that, 
while the BCS1 and RLN gaps are smaller than the BCS ones at a given G, especially for 
the BCS1 gap at G — G^. CS1 , the increases of the gap offered by the LN1 and RLN compared 
to the LN value are negligible at all G. 



B. Ground-state energy 

Shown in Fig. [2] are the results for the ground-state energies obtained within the BCS, 
LN, SCRPA, QRPA, LNQRPA, SCQRPA, and LNSCQRPA in comparison with the exact 
one for N = 10. The exact result is obtained by directly diagonalizing the Hamiltonian 
in the Fock space 21] . It is seen that the BCS strongly overestimates the exact solution. 
The LN result comes much closer to the exact one even in the vicinity of the BCS (QRPA) 
critical point, while the QRPA (RPA) result agrees well with the exact solution only at 
G ^> G„ cs (G <C G^ r cs ). The improvement given by the SCRPA is significant as its result 
nearly coincides with the exact one in the weak coupling region. However the convergence 
of the SCRPA solution is getting poor in the strong coupling region. As a result, only 
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FIG. 2: (Color on line) Ground state energies as functions of G for N = 10. The exact result is 
represented by the thin solid line in both panels (a) and (b). In panel (a), the dotted line denotes 
the BCS result, the thin dashed line stands for the LN result, the dash-dotted line shows the pp 
RPA result at G < G* cs , and the QRPA one at G > G® cs , while the dash-double-dotted line 
depicts the LNQRPA result. Predictions by self-consistent approaches are plotted in panel (b), 
where the thick dashed line denotes the SCRPA result, while the SCQRPA and LNSCQRPA are 
shown by the thick solid and double-dash-dotted lines, respectively. 

the values up to G < 0.46 MeV are accessible. The SCQRPA is much better than the 
QRPA as it fits well the exact ground-state energy at G > G^ CS1 . The LNQRPA strongly 
underestimates the exact solution while the LNSCQRPA, which includes the effects due to 
the screening factors in combination with PNP, significantly improves the overall fit. From 
this analysis, we can say that, among all the approximations undergoing the test to describe 
simultaneously the ground and excited states, the SCRPA, SCQRPA, and LNSQRPA can 
be selected as those which fit best the exact ground-state energy. The LN method based 
on the BCS (thin dashed line) also fits quite well the exact one at all G but it does not 
allow to describe the excited states as the approaches based on the QRPA do. Although the 
fit offered by the LNSCQRPA in the vicinity of the critical point is somewhat poorer than 
those given by the SCRPA and the SCQRPA, its advantage is that it does not suffer any 
phase-transition point due to the violation of particle number as well as the Pauli principle. 
The corrections due to ground-state correlations can also be clearly seen by examining 
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TABLE I: The energy difference AE = E g , s ,(G) - E g . B .(0) at various G (in MeV) as predicted by 
the QRPA, SCQRPA, LNQRPA, LNSCQRPA, and exact solutions for N = 10. 



G 


QRPA 


SCQRPA 


LNQRPA 


LNSCQRPA 


Exact 


0.10 






-0.05 


-0.06 


-0.04 


0.20 






-0.24 


-0.28 


-0.17 


0.30 






-0.63 


-0.69 


-0.44 


0.35 


-0.93 




-0.91 


-0.94 


-0.64 


0.40 


-1.00 




-1.26 


-1.21 


-0.90 


0.47 


-1.38 


-1.44 


-1.86 


-1.66 


-1.36 


0.50 


-1.60 


-1.66 


-2.16 


-1.88 


-1.60 


0.60 


-2.53 


-2.58 


-3.34 


-2.80 


-2.56 


0.70 


-3.70 


-3.75 


-4.76 


-3.96 


-3.76 


0.80 


-5.09 


-5.13 


-6.39 


-5.33 


-5.17 


o.yo 


-6.65 


-6.68 


-8.19 


-6.87 


-0.75 


1.00 


-8.34 


-8.38 


-10.13 


-8.56 


-8.46 


1.10 


-10.15 


-10.18 


-12.19 


-10.37 


-10.29 


1.20 


-12.05 


-12.08 


-14.33 


-12.27 


-12.22 


1.30 


-14.03 


-14.06 


-16.55 


-14.25 


-14.22 


1.40 


-16.06 


-16.10 


-18.84 


-16.30 


-16.28 



the energy difference 

AE = E g . s .(G) - E g . s .(0) (63) 

between the ground-state energies defined at finite and zero G 2 . The values of this energy 
difference as predicted by the QRPA, SCQRPA, LNQRPA, and LNSCQRPA for the system 
with N = 10 at various G are compared with the exact ones in Table [B It is seen from 

2 Within the RPA and SCRPA, where the mean field is the HF one, AE coincides with the correlation energy 
E COIT = E s , s . — -Ehf because Eg C .s act \o) = £"hf 7 (f^ F = 0, //f F = 1). Within the quasiparticle formalism, 
however, B corr is defined as the difference between the QRPA (LNQRPA, SCQRPA, LNSQRPA) ground- 
state energy and that given within the BCS (LN, LN1) method. This -Ecorr is quite different from AE in 
the strong-coupling regime because of the large pairing gap. Therefore we find more appropriate in the 
quasiparticle representation to compare the approximated and exact energies AE (1551) rather than E COTI . 



16 



TABLE II: Relative errors SE^ and 5E^ from Eq. flSU) at various G as predicted by the QRPA, 
SCQRPA, LNQRPA, and LNSCQRPA for N = 10. 





6E& (%) 




5E^ 


(%) 




vj iivie v i 


QRPA SCQRPA LNQRPA LNSCQRPA 


QRPA SCQRPA LNQRPA LNSCQRPA 


n 1 n 

U. 1U 




25.00 


50.00 






0.04 


0.08 


u.zu 




41.18 


64.71 






0.28 


0.44 


U.OU 




43.18 


56.82 






0.75 


0.98 


n ^ 

U.OtJ 


43.51 


42.19 


46.88 


1.13 




1.05 


1.17 


U.4iU 


11.11 


40.00 


34.44 


0.39 




1.39 


1.20 


n 47 


1.47 


5.88 36.76 


22.06 


0.08 


0.30 


1.90 


1.14 


U.OU 


0.00 


3.75 35.00 


17.50 


0.00 


0.23 


2.11 


1.05 


n fin 
u.uu 


1.17 


0.78 30.47 


9.37 


0.11 


0.07 


2.83 


0.87 


u. / u 


i fin 


n 07 ofi fin 


o.oz 


U.zl 


n hq 


Q A Q 

o.4o 


n 7n 
U. 


U.OU 


1.55 


0.77 23.60 


3.09 


0.27 


0.13 


4.04 


0.53 


0.90 


1.48 


1.04 21.33 


1.78 


0.32 


0.22 


4.54 


0.38 


1.00 


1.42 


0.95 19.74 


1.18 


0.36 


0.24 


4.99 


0.30 


1.10 


1.36 


1.07 18.46 


0.78 


0.40 


0.31 


5.38 


0.23 


1.20 


1.39 


1.15 17.27 


0.41 


0.46 


0.38 


5.67 


0.13 


1.30 


1.34 


1.13 16.39 


0.21 


0.48 


0.41 


5.94 


0.08 


1.40 


1.35 


1.11 15.72 


0.12 


0.53 


0.44 


6.20 


0.05 



this table that, while in the weak coupling regime (Gjf r cs < G < 0.8 MeV) the QRPA and 
SCQRPA predictions for this energy difference are closer to the exact result, at high G the 
SCQRPA and LNSCQRPA are the ones that offer the better fits for this quantity. The 
LNQRPA, on the contrary, offers a quite poor fit for AE to the exact result. 
A more quantitative calibrations can be seen by analyzing the relative errors 

A p(approx) A p(cxact) 77>(approx) p(exact) 

<*# (a) = — . / , and 5E^ = , (64) 

A TP (exact) TP (exact) v ' 

which are shown in Table [Til Because A£^ exact - ) are quite small at small G, the relative 
errors 5E^ are quite large in the weak-coupling region. In this respect the relative error 
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5E^ turns out to be a better calibration. While 5E^ decreases as G increases for all 
approximations with the LNSCQRPA having the smallest relative errors at large G, the 
behavior of SE^ on G is somewhat different depending on the approximation. A decrease 
of this quantity is seen within the QRPA and SCQRPA with increasing G up to G = 0.7 
MeV, and an increase with G takes place at large G. For the LNSCQRPA, the relative 
error SE^ increases first with G up to G = 0.4 MeV, then decreases at larger G. Within 
LNQRPA one sees a steady increase of 5E^ with G to reach a value as large as 6.2 % at 
G = 1.4 MeV. 

The quantities that are directly defined by the differences of ground-state energies are 
the chemical potentials A ± and A, namely 

A+ = i [E g . s .(N + 2) - E S . S .(N)} , A- = ~ [E g . s .(N) - E g . s .(N - 2)] , A = ^(A+ + A~) . 

(65) 

The exact values of the chemical potentials A and A 1 * 1 are shown in Fig. [3] in comparison with 
the predictions within quasiparticle presentations for N = 10. It is seen from this figure that 
the SCRPA and SCQRPA [Fig. EJ(d) -13(f)] offer the best fit to the exact results except that 
the SCRPA poorly converges at G > 0.4 MeV, while SCQRPA stops at G = G£ CS1 . The 
RPA and QRPA also describe very well the exact results, except the values in the critical 
region, where the RPA and QRPA diverge. The LNSCQRPA predictions for the chemical 
potentials show smooth functions at all G, which fit well the exact results, including the 
region around G cr , where they slightly underestimates the exact ones. 



C. Energies of excited state 

As has been discussed in Refs. Q, the first solution u x of the QRPA or SCQRPA 
equations is the energy of spurious mode, which is well separated from the physical solutions 
uo v with v > 2. The first excited state energy is therefore given by u>2- Figure H] shows the 
exact eigenvalues for the excited states. As has also been demonstrated in Ref. 22j, this 
figure shows that the coupling in the small-G region causes only small perturbations in the 
single-particle levels. With increasing G the system goes to the crossover regime, where 
level splitting and crossing are seen, releasing the levels' degeneracy. In the strong coupling 
regime the levels coalesce into narrow well-separated bands. The approaches based on the 
QRPA with PNP within the LN method also splits the levels but the nature of the splitting 

18 



-5.0 

_ -5.2 
> 

3 -5.4 



-5.6 
-5.8 



-5.6 
© -5.8 
— -6.0 
^ -6.2 

-6.4 



> 

J 4.8 



I -5.0 
-5.2 
-5.4 



1 1 1 1 1 1 1 

' n = io : 

- (a) Y 

— i i — i — i — i — i — i — 


I 1 1 1 1 1 1 

■ (d) V 


- (b) \ 


— i — i — . — i — . — i — i — 
■ (e) \ 


— i — i — i — i — i — i — i — 
- (c) \ 

I.I.I. 1 


— i — i — i — i — i — i — . — 

■ (f) \ 



0.2 0.4 0.6 
G (MeV) 



0.2 0.4 0.6 
G (MeV) 



0.8 



FIG. 3: (Color online) Chemical potentials A and X^ 1 as functions of G for N = 10 as predicted by 
the exact solutions, RPA, QRPA, SCRPA, SCQRPA, and LNSCQRPA. Notations are as in Fig. 
2. 

comes from the two components within the QRPA operator fHUj) . which correspond to the 
addition and removal modes, respectively, in the RPA limit. When the pairing gap A is 
finite, it is not possible to consider the QRPA excitations as purely addition or removal 
modes, but only as those with some components having the dominating property inherent 
to one of these modes. The QRPA eigenvalues also have two branches with positive u u and 
negative — lo v energies. However, unlike the pp RPA, where the negative eigenvalues in the 
equations for addition modes are also physical as they are the energies of the removal modes 
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FIG. 4: Exact energies E^ x = £^ X (N) — £q x (N) obtained within the Richardson model for excited 
states v relative to the exact ground-state level £q x as functions of G for N=10. 

taken with the minus sign and vice versa, within the QRPA only the positive energies uj v are 
physical, and they are compared with the exact ones, E^ x = £ U (N) — £ (N), in the present 
paper. 

As an example to illustrate this level-splitting pattern, we show in Fig. [5] the exact energy 
Ef 1 = £i(N) — £o(N) of the lowest excited state [y = 1) with respect to the exact ground 
state [y = 0) in the system with N = 10 particles in comparison with the predictions within 
the QRPA, LNQRPA, SCQRPA, and LNSCQRPA 3 . As the exact energy represents 
the energy of the lowest pair-vibration state, it is compared with the energies u\ of the 
lowest excited state obtained within QRPA, LNQRPA, SCQRPA and LNSCQRPA, which 
are built on the pairing condensate (quasiparticle vacuum). The splitting is clearly seen 
from Fig. [5] within the LN method, namely the LNQRPA and LNSCQRPA. One can see 
that, within the LN(SC)QRPA, each single level at G = splits into two components in the 

,i n 41 LNQRPA , LNQRPA LNSCQRPA , LNSCQRPA rp , i 

small-C region, e.g. the pair u 2 and u 3 orw 2 and uj 3 ^ . lo look 

inside the source of the splitting, we rewrite the QRPA operator ( HO"]) into the components 



3 For the two- level case E^ x corresponds to the solid line in the upper panel of Figs. 1, 3 - 5 in Ref. [9| or 
Figs. 1 - 3 in Ref. ',17\ for N = 4, 8, and 12). 
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FIG. 5: (Color on line) The energies of the first excited state as functions of G at A=10. The 
results refer to the exact solution, E^ x (solid line), the QRPA solution, u^ RPA (dash-dotted line), 
the SCQRPA solution, u>2 C( °^ RPA (thick solid line), the LNQRPA solutions, u>2 N( ^ RPA (thin dash 
- double-dotted line) and u;3 N ^ RPA (thick dash - double-dotted line), as well as the LNSCQRPA 
solutions, ^LNSCQRPA double-dash - dotted line) and ^^nscqrpa ^ick double-dash - dotted 
line). 

with dominating contributions of addition- and removal-mode patterns as follows: 

Ql = (Ql) {A) + (Ql) {R) , 

(QiY A) = E x p a I - E > (Qt ) (R) = E x * A i - E Y v A v , (66) 

p h hp 

where the indices j run over all the levels, from which those located below (above) the 
chemical potential are formally labelled with h (p) indices. It is not difficult to see that, in 
the RPA limit (or zero-pairing limit), (Qt)^ is transformed into operator A\, that generates 
the addition modes, while (Ql)^ becomes B) v that generates the removal modes (in the 
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FIG. 6: The energies of the first excited state in different schemes as functions of G for N = 10. 
The thin and thick dash - double-dotted lines denote the second and third LNQRPA solutions, 
while the thin and thick dotted lines stand for the absolute values of the corresponding solutions 
within the LNQRPA1 scheme. 

standard notations for addition and removal operators from Refs. jj, 5, 6]). Using this formal 
expression ( |66l) . we derived the QRPA equations for the excitations generated by operators 
(Qt) (A) and {Ql) (K \ separately. The energies of the corresponding first excited states from 
the resulting sets of equations were calculated by using the LN method. We call this scheme 
as LNQRPA1. The set of equations for the modes generated by operator (Ql)^ gives a 
negative u;2 N< ^ RPA1 and positive u>3 N( ^ RPA1 , which means that they correspond to the energies 
of the removal and addition modes, respectively. The absolute values of these energies are 
shown in Fig. [6] along with u\ g^ RPA . It is seen from this figure that in the weak-coupling 
region the higher-lying levels ClJ ,^ jN Q R ' PA anc i ^lnqrpai near iy coincide, while the lower-lying 
one, d;2 N( ^ RPA , is almost the same as |a>2 N( ^ RPA1 |. From here, we can identify u; LN( ^ RPA and 
^lnqrpa ag i eve i g wnere the addition and removal modes dominate, respectively. As the 
interaction G increases, the occupation probabilities of the levels below and above the Fermi 
level become comparable so it becomes more and more difficult to separate the patterns 
belonging to addition and removal modes in the QRPA excitations. 

From this analysis and Fig. [5], it becomes clear that, in the weak coupling region, the 
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level o;3 N( ^ RPA , which is generated mainly by the addition mode, fits well the exact result, 
while the agreement between the exact energy and cu^ RPA as well as co»2 C< ^ RPA is good only 
in the strong coupling region. At large values of G, predictions by all approximations and 
the exact solution coalesce into one band, whose width vanishes in the limit G — > oo. 

The energies of the ground state and the first excited state obtained for TV = 4, 6, 8 are 
depicted in Fig. [7J The figure shows that increasing TV worsens the agreement of the results 
obtained within the LNQRPA and LNSCQPPA with the exact ones for both the ground 
state and the first excited state, while the QRPA and SCQRPA results become closer to the 
exact ones at G 3> G CT . At small TV (TV = 4), the solution u;3 N< ^ RPA seems to fit best the 
exact result for all values of G. 

The pair-vibration excitation energy Ef x is usually larger than the energy of the lowest 
state with one broken pair. The latter is described within the pp RPA as the energy of 
the lowest addition mode in the laboratory reference frame fixed to the ground state of 
TV-particle system . It is worthwhile to compare the predictions for the excited-state 

energies obtained within the quasiparticle approaches developed in the present paper with 
pp RPA and SCRPA predictions by transforming the latter into the intrinsic reference frame 
of the system with TV + 2 particles. This is done as follows. From the (SC)RPA energy of the 
ground-state level u;q SC ' )RPA = £q SC ^ RPA (TV + 2) — £q SC ^ rpa (TV), and that of the first excited 
state w f °) RPA = £f C > RPA (TV + 2) - £<f C)RPA (TV) 4 it follows that 

A ^(SC)RPA s ^(SC)RPA _ ^SQRPA = ^(SQRPA^ + 2) _ ^(SC)RPA (iv + 2) ? ^ 

This energy Aa/ SC ) RPA is shown in Fig. [8] as a function of G along with the corresponding 
LNQRPA, LNSCQRPA, and exact energies for several values of TV. This figure clearly shows 
that the LNQRPA and LNSCQRPA are superior to the pp RPA and SCRPA as they offer 
an overall prediction closer to the exact result for all G and TV. They neither collapse at 
a G cr as in the case with the pp RPA nor have a poor convergence as the SCRPA does at 
G > G cr . 



4 The energies lJq SC ^ RPA and cj 1 SC ^ RFA correspond to energies E\ and E 2 shown in Figs. 3 and 4 in Ref. 
respectively. 
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G (MeV) G (MeV) 

FIG. 7: (Color on line) Energies of ground state (left panels) (notations as in Fig. [2|) and first 
excited state (right panels) (notations as in Fig. [5]) for several values of iV indicated on the panels 
as functions of G. 

IV. CONCLUSIONS 

This work proposes a self-consistent version of the QRPA in combination with particle- 
number projection within the Lipkin-Nogami method as an approach that works at any val- 



24 



1 




(c) N = 6 



1.5 



1.0 



(d) N = 4 



0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 

G (MeV) G (MeV) 



FIG. 8: (Color on line) Energy Au/ SC ) RPA (jSTJ) obtained within the pp RPA (dash-dotted line) 
and SCRPA (thick solid line) as a function of G for several values of N in comparison with the 
energy u;g N( ^ RPA (dash - double-dotted line), ^lnscqrpa (double-dash - dotted line), and the 
exact energy Ef x (thin solid line), which are the same as those in Fig. [7] (d) - (f) for N = 4, 6, 
and 8. 

ues of the pairing-interaction parameter G without suffering a phase-transition-like collapse 
(or poor convergence) due to the violation of Pauli principle as well as of the integral of mo- 
tion such as the particle number. The self-consistency is maintained within a set of coupled 
equations for the pairing gap, QRPA amplitudes, and energies by means of the screening 
factors, which are the expectation values of the products of quasiparticle-pair operators, 
and the ground-state correlation factor, which is a function of the QRPA backward-going 
amplitudes. 

The proposed approach is tested in a multi-level exactly solvable model, namely the 
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Richardson model for pairing. The energies of the ground and first-excited states are calcu- 
lated within several approximations such as the BCS, RBCS, BCS1, LN, RLN, LN1, QRPA, 
SCQRPA, LNQRPA and LNSCQRPA. The obtained results for the ground-state energy 
show that the use of the LN method that includes the SCQRPA correlations not only allows 
us to avoid the collapse of the BCS as well as the QRPA but also fits well the exact result. 
For the energy of the first excited state, the LNQRPA and LNSCQRPA results offer the 
best fits to the exact solutions in the weak coupling region with large particle numbers, 
while the QRPA and SCQRPA reproduce well the exact one in the strong coupling region. 
In the limit of very large G all the approximations predict nearly the same value as that of 
the exact one. As the number of particles decreases, it becomes sufficiently well to use the 
predictions given by the LNQRPA and LNSCQRPA for energies of both the ground state 
and first-excited state to fit the exact results. 

We believe that the approach proposed in this work can be useful in the applications to 
light and unstable nuclei, where the validity of the QBA and that of the conventional BCS 
are in question. Such applications are the goal for forthcoming studies. 
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APPENDIX A: ACCURACY OF APPROXIMATION (f33|) 

Let us analyze the accuracy of the assumption fl33|) used in the numerical solutions of the 
BCS1, LN1, and SCQRPA equations in the present paper. 

Shown in the 2nd and 5th columns of Table II I II are the values of the pairing gaps A and 
A obtained under the approximation (1331) within the BCS1 and LN1 method, respectively. 
They are compared with the average gaps A (3rd column) and A (6th column), which are 
the values obtained by averaging the level-dependent BCS1 gap Aj and LN1 gap Aj over 
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TABLE III: BCS1 and LN1 pairing gaps (in MeV) at various values of G (in MeV) (see text). 





BCS1 


LN1 


G 


A 


A 


% (%) 




A 




A 


f (%) 





01 











0015 





0015 


0.0000 





10 











0606 





0607 


0.1647 





20 











2279 





2289 


0.4369 





30 











5278 





5321 


0.8081 





40 











9579 





9660 


0.8385 





47 


0.8224 


0.8357 


1.5915 


1 


3139 


1 


3233 


0.7103 





50 


1.0694 


1.0829 


1.2467 


1 


4742 


1 


4839 


0.6537 





60 


1.7219 


1.7351 


0.7608 


2 


0261 


2 


0360 


0.4862 





70 


2.3314 


2.3436 


0.5206 


2 


5896 


2 


5993 


0.3617 





80 


2.9279 


2.9391 


0.3811 


3 


1541 


3 


1633 


0.2908 





90 


3.5132 


3.5235 


0.2923 


3 


7148 


3 


7234 


0.2310 


1 


00 


4.0882 


4.0977 


0.2318 


4 


2701 


4 


2783 


0.1917 


1 


10 


4.6539 


4.6629 


0.1930 


4 


8197 


4 


8277 


0.1657 


1 


20 


5.2118 


5.2203 


0.1628 


5 


3641 


5 


3718 


0.1433 


1 


30 


5.7628 


5.7710 


0.1421 


5 


9037 


5 


9113 


0.1286 


1 


40 


6.3079 


6.3160 


0.1282 


6 


4390 


6 


4466 


0.1179 



all the levels, namely A = ^ . Aj/N and A = £\ Aj/N. The secon d term at the rhs of Eq. 
(I31~j) . which contains 8J\fj as evaluated by the approximation (1321 . is taken into account in 
calculating Aj and Aj within the perturbation theory, i.e. with rij being evaluated within 
SCQRPA and LNSCQRPA (where this term is neglected). Except for the two values at 
G = Gjf r csl = 0.47 MeV and G = 0.5 MeV within the BCS1, we see that the values of the 
relative errors 8 A/ A = (A — A) /A and 5 A/ A = (A — A)/A are all smaller than 1 %, and 
decrease with increasing G. 

Shown in Table HVl are the values of the ratio (5A/}) 2 / (T>j) from Eqs. (|31~1) and fl32|) corre- 
sponding to the five lowest levels for N = 10 at various G obtained within the LNSCQRPA. 
The largest value of this ratio is observed at the level with j = 5, the closest one to the Fermi 
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TABLE IV: The ratio (5J\fj) 2 / '(Pj) from Eqs. (|3ip and (|32p corresponding to the 5 lowest levels 
j = 1, . . . , 5, and the energies 1^3 (in MeV) of the first excited state described in the text for 
N = 10 at different values of G (in MeV) within the LNSCQRPA. The energy oo^{a) is obtained 
including the last term at the rhs of Eq. (|3ip . while 0^3(6) is calculated using the approximation 



G 


3 = 1 


3 = 2 


3 = 3 


j = 4 


3 = 5 


w 3 (a) 


u 3 {b) 


0.01 


0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


2.0001 


2.0001 


0.2 


0.0009 


0.0012 


0.0017 


0.0027 


0.0046 


2.0697 


2.0711 


0.4 


0.0023 


0.0030 


0.0040 


0.0055 


0.0082 


2.6701 


2.6742 


0.6 


0.0019 


0.0023 


0.0027 


0.0032 


0.0054 


4.2040 


4.2067 


0.8 


0.0013 


0.0015 


0.0016 


0.0021 


0.0033 


6.1514 


6.1531 


1.0 


0.0009 


0.0010 


0.0011 


0.0015 


0.0022 


8.1798 


8.1812 


1.2 


0.0006 


0.0007 


0.0009 


0.0012 


0.0017 


10.211 


10.212 


1.4 


0.0005 


0.0006 


0.0008 


0.0011 


0.0014 


12.229 


12.230 



level, at G = 0.4 MeV (close to G^ GSl ). But it amounts to only 0.0082, which is a clear 
evidence that this ratio is indeed negligible. The last two columns of this table display the 
energies uJs(a), obtained within the LNSCQRPA including the last term at the rhs of Eq. 
(I31~]) . and 003(b), which the LNSCQRPA predicts within the approximation (|33|) . Although a 
systematic 0^3(0) > 003(b) is observed, the largest difference, also seen at G = 0.4 MeV, does 
not exceed 0.15 %. These results guarantee the high accuracy of the approximation (|33l) . 
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